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We study the effects of generic short-ranged interactions on a system of 2D Dirac fermions subject 
to a special kind of static disorder, often referred to as "chiral." The non-interacting system is a 
member of the disorder class BDl [M. R. Zirnbauer, J. Math. Phys. 37, 4986 (1996)]. It emerges, 
for example, as a low-energy description of a time-reversal invariant tight-binding model of spinless 
fermions on a honeycomb lattice, subject to random hopping, and possessing particle-hole symmetry. 
It is known that, in the absence of interactions, this disordered system is special in that it does 
not localize in 2D, but possesses extended states and a finite conductivity at zero energy, as well 
as a strongly divergent low-energy density of states. In the context of the hopping model, the 
short-range interactions that we consider are particle-hole symmetric density-density interactions. 
Using a perturbative one-loop renormalization group analysis, we show that the same mechanism 
responsible for the divergence of the density of states in the non-interacting system leads to an 
instability, in which the interactions are driven strongly relevant by the disorder. This result should 
be contrasted with the limit of clean Dirac fermions in 2D, which is stable against the inclusion 
of weak short-ranged interactions. Our work suggests a novel mechanism wherein a clean system, 
initially insensitive to interaction effects, can be made unstable to interactions upon the inclusion 
of weak static disorder. We dub this mechanism a "disorder-driven Mott transition." Our result for 
2D fermions also contrasts sharply with known results in ID, where a similar delocalized phase has 
been shown to be robust against the inclusion of weak interaction effects. 
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I. INTRODUCTION 

The discovery^ of novel symmetry classes of disor- 
dered quantum systems has reawakened interest in ( "An- 
derson") localization physics. These classes describe the 
physics of non-interacting quantum particles subject to 
static disorder potentials on scales much longer than the 
mean free path. Model realizations include descriptions 
of superconducting quasiparticles^iiiSiSiLSiS, as well as 
models of electrons with random hopping on bipartite 
lattices ("sublattice models" ) 10 i 11 i 12 i 13 i 14 i 15 . In contrast 
to the conventional disordered metals, characterized by 
the standard and much studied ( "unitary" , "orthogonal," 
and "symplectic" ) Wigner-Dyson classes, models in these 
new classes often exhibit surprising low-dimensional be- 
havior. Some such disordered systems possess delocalized 
states in one and two dimensions, and these states can 
persist even in the presence of strong randomness. (See 
e.g. Refs. [l0llTlTr| .) 

While much progress has been made in understanding 
the disorder physics in such systems of non-interacting 
particles, the role of interparticle interactions is often 
poorly understood. This issue is crucial to the search 
for experimental realizations. In systems belonging 
to the conventional (metallic) Wigner-Dyson classes, it 
is knownii that interaction effects can drastically al- 
ter the localization transitions. For example, in all of 
these metallic universality classes, the incorporation of 
long-range Coulomb interactions destabilizes the non- 
interacting Anderson metal-insulator transition; in its 
place, a new, interacting "Anderson-Mott" transition is 
typically found in D = 2 + e dimensions, 17 character- 
ized by distinct critical exponents. Renewed interest in 



these issues was recently generated through discussions 
concerning a metal-insulator transition in 2Di& 

In this paper, we examine the effects of generic short- 
ranged interactions upon a 2D example of the so-called 
"sublattice" (or: "chiral") class denoted by BDl (in the 
classification scheme of Ref. [lj). As a particular repre- 
sentative of this class we study a time-reversal invariant 
(TRI) model of spinless fermions with random nearest- 
neighbor hopping on the honeycomb lattice. The random 
single-particle Hamiltonian h is special because the hop- 
ping occurs only between the two triangular sublattices 
of the honeycomb lattice [Eq.@ below]. This translates 
into an additional, discrete, so-called "chiral" symmetry, 
present in every realization of disorder 

-a 3 ho- 3 = h, (1) 

where <r 3 is a Pauli matrix diagonal in sublattice space. 
The chiral symmetry can be thought of as a combination 
of time-reversal and particle-hole symmetries. The hon- 
eycomb model admits a low-energy description in terms 
of Dirac fermions subject to two types of long-wavelength 
random potentials: a pair of random masses, and a ran- 
dom vector potential. The resulting disordered Dirac 
effective field theory also applies to the random bond 
7r-flux model on the square lattice ^Li2iS2i The theory is 
describe d 11 ' 13 by two independent parameters, which we 
will call g and gA, characterizing the sample-to-sample 
fluctuations of the zero mean Gaussian random mass 
and random vector potential disorders, respectively (see 
Sec. Ill Al for a review). 

A key theoretical motivation for studying the disor- 
dered Dirac model with chiral symmetry is that it consti- 
tutes one of the simplest known arenas for investigating 



several non-trivial aspects of disorder-dominated quan- 
tum criticality. While some features of the physics of 
the sublattice model are certainly tied to the special chi- 
ral symmetry, others, particularly the multifractal na- 
ture of the extended wavefunctionsiiSSaSiiSSiSi 2 ^ 5 . are 
believed to be more general features of (Ander son) lo- 
calization critical points (see e.g. Refs. |2ll22|| ). Chi- 
ral disorder models exhibit the most interesting physics 
in one and two dimensions, where these systems have 
proven amenable to a variety of analytical techniques, 
and several exact and/or non-perturbative results are 
known 13 i 14 i 15 i 16 i 26 i 27 . 

Models in the chiral symmetry classes were originally 
studied analytically by Gade and Wegneui ", who formu- 
lated a non-linear sigma model (NLcrM) appropriate to 
a system of non-interacting fermions subject to chiral 
disorder! 41 . Gade and Wegner studied cases both with 
and without time-reversal invariance in the weak cou- 
pling limit and found essentially identical physics. Just 
as the Dirac theory mentioned above, the NLerM also 
contains two parameters: the first is the usual dimen- 
sionless DC resistance, which is associated with the ran- 
dom mass coupling g in our Dirac theory. The param- 
eter gA, on the other hand, is associated with an addi- 
tional "U(l)-term" , which is present only in the chiral 
class NLcrMs. A perturbative one-loop renormalization 
group (RG) analysis performed in Ref. ^3 yielded that, 
in 2D, the DC resistance does not renormalize, while the 
strength of the "U(l)-term" is driven to strong coupling. 
Numerical work on 2D chiral models includes that listed 
e.g. in 

The TRI chiral Dirac theory studied in the present pa- 
per was introduced by Hatsugai, Wen, and Kohmoto in 
Ref. Ull where results of a one-loop RG analysis, adapted 
from Ref. 0) were presented. Using techniques from con- 
formal field theory, Guruswamy, LeClair, and Ludwig 13 
later performed a non-perturbative analysis (in g and 
gA)- In particular, closely parallel to the results for the 
NLerM discussed above, the strength g of the random 
mass is an exactly marginal perturbation (in the RG 
sense) to the clean Dirac theory, while the strength gA of 
the random vector potential is driven to strong coupling 
by the presence of a non-zero g. The pure marginality of 
g indicates the existence of a critical, delocalized phase 
at zero energy, made possible by the special chiral sym- 
metry. The running gA parameter leads to the prediction 
of a strongly divergent low-energy density of states 

v(uj) ~ — exp (— c |lnw| Q ) , (2) 

where the scale-independent constant c is a function of g 
and the exponent a = 1/2. 

More recently, Motrunich, Damle, and Huse 1 ^ used a 
strong randomness RG picture to argue that, at asymp- 
totically low energies, or equivalently, large values of the 
parameter gA, the density of states of the chiral Dirac 
model should retain the form as in Eq. @, but with a 
different exponent a — 2/3. Mudry, Ryu, and Furusak: 15 
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FIG. 1: The honeycomb lattice. 



observed that the freezing transition 2 * 2 ^, which was 
known to occur for vanishing g (the abelian gauge poten- 
tial mode l 20 ! 35 ) when gA exceeds some critical value g\, 
affects the dynamic critical exponent z, and thus the low- 
energy density of states. They noted that the previous 
analytical RG treatments did not take into account an 
infinite set of operators associated with moments of the 
local density of states (LDOS). These operators carry in- 
creasingly negative (i.e. increasingly relevant) scaling di- 
mensions, whose appearance is taken to signal the broad- 
ening of the probability distribution function for the 
LDOS. Using a functional renormalization group treat- 
ment pioneered by Carpcntier and Le Doussal^S to or- 
ganize the diverging moments of the LDOS, the authors 
of Ref. 10 derived the dynamic critical exponent analyt- 
ically and thereby confirmed the asymptotic behavior of 
the global density of states with a = 2/3 of Ref. IbH 

A second motivating factor for our work is that the 
ID version of our 2D sublattice model studied in this 
paper is known to be stable to the inclusion of weak 
short-ranged interactions (which are irrelevant in the RG 
sense) 16 - 26 . The former model also possesses a strongly- 
divergent low-energy density of states, with a critically 
delocalized state at zero energyiSJ. The robust nature of 
the ID critical phase to the inclusion of interaction effects 
is somewhat counterintuitive, given the picture one has 
of a large number of states crammed into a very narrow 
window of energy. In this paper, we will show that this 
result does not generalize to 2D. 

In the present paper, we add to the continuum Dirac 
description of the disordered honeycomb hopping model 
generic short-ranged interactions that preserve the chiral 
symmetry, defined for the Dirac model below in Eq. JSJl, 
and in the paragraph following that equation. We use the 
Schwinger-Keldysh method 37 to perform the ensemble 
average over realizations of the disorder. We then employ 
a one-loop renormalization group treatment to investi- 
gate the stability of the disordered, critically delocalized 
phase to interaction effects. We find that the same mech- 
anism responsible for driving the divergence of the low- 
energy density of states also feeds into same-sublattice 
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(e.g. next-nearest-neighbor) interactions, producing an 
instability where the interaction strength grows to large 
values. We show that this instability always occurs at 
low enough energies, regardless of whether or not the 
(non-interacting) disordered system reaches the regime 
beyond the freezing transition discussed above. We spec- 
ulate that this instability corresponds to the onset of 
some kind of Mott insulating phase, dominated by in- 
teractions rather than disorder. 

Our results are both expected and surprising: ex- 
pected, because of the strongly divergent density of states 
and the potentially fragile nature of the low-dimensional, 
critical disorder-only physics (i.e. in the absence of inter- 
actions); surprising, first because, as was already men- 
tioned, the opposite result obtains in ID, and second be- 
cause the clean Dirac fixed point is robust against all 
short-ranged interaction effects. If the observed insta- 
bility ultimately drives the disordered chiral symmetric 
model into an interaction-dominated, homogeneous Mott 
insulating phase, then we have discovered a somewhat cu- 
rious route: a clean system, initially insensitive to inter- 
action effects, can be made unstable to interactions upon 
the inclusion of (a special kind of) weak static disorder. 
One might dub this a "disorder-driven Mott transition." 

The organization of this paper is as follows. In Sec. ITT1 
we define the honeycomb lattice model and summarize 
the disordered low-energy effective Dirac theory. We then 
include short-ranged interactions, and define a set of cou- 
pling constants to encode the model physics. In Sec. IIIII 
we set up a one-loop renormalization group treatment of 
the full, disordered and interacting field theory. We re- 
derive the weak-coupling RG flow equations of the non- 
interacting, disorder-only modeliAai&M as a check of our 
methodology. We then compute the scaling dimensions 
of the various interaction operators; these scaling dimen- 
sions are new. In Sec. IIVI we solve the combined RG 
flow equations for the disorder and interaction couplings, 
and we interpret the results. The reader less interested in 
calculational details may skip Sec. 11111 and proceed from 
the end of Sec. HJ directly to Sec.HVl 



II. MODEL AND SETUP 

A. Lattice model and continuum limit 

Consider a system of spinless fermions at half-filling 
on the honeycomb lattice with nearest-neighbor hopping. 
The lattice is depicted in Fig. ^ We define fermion anni- 
hilation operators ca% and CBj on the "A" and "B" trian- 
gular sublattices, respectively; then the random hopping 
Hamiltonian may be written as 



H 



E4 



Uj cb 3 + h.c. 



(3) 



where the sum runs over all nearest-neighbor bonds on 
the honeycomb lattice. We take the hopping amplitudes 



tij to be purely real. In the clean limit, tjj = t; then 
the Hamiltonian in Eq. possesses two energy bands 
of Bloch states indexed over the hexagonal Brillouin zone 
(Fig. |2J), and, as in graphene, the low energy states reside 
near the corners of the zone where the two bands meet 
at isolated Fermi points. 

We truncate the Bloch spectrum of the clean system 
to low-energy modes with wavenumbers k — (k x ,k y ) 
living within a momentum cutoff A of the zone cor- 
ners. Alternate corners are equivalent by reciprocal lat- 
tice translations, so we retain modes near two oppo- 
site corners, which we index by the Fermi wavevectors 
±k F = (±k F ,0), shown in Fig. An effective low- 
energy Dirac description emerges if we construct the 4- 
component spinor 



c A (k + k F ) 
c B (k + k F ) 
CA(k - k F ) 
c B (k - k F ) 



(4) 



with |fc| < A. Linearizing the energy bands near the 
Fermi points and incorporating the hopping disorder as 
a perturbation, we find the non-interacting Hamiltonian 



H 



D 



d 2 rip\r) [v F (id^ 
d 2 rijj\r) [m % l3 l ] ip(r 



A M j9°)]V(r) 



(5) 



where v F is the Fermi velocity, and ij)(r) is the long- 
wavelength Fourier transform of J3J. The disorder 
emerges in the continuum theory ||SJ| as a random vector 
potential A^ix), fj, £ and a pair of random Dirac 

masses m l (x), i € {1,2}. In Eq. JSJ, summation on re- 
peated indices is implied and will be assumed throughout 
the rest of the paper. 

The symbols (a x ,a y ) and {/3°, /3 1 , f3 2 } in Eq. JSJ) de- 
note various 4x4 matrices. Explicit representations can 
be derived from the honeycomb lattice model. We in- 
troduce two sets of Pauli matrices^ the matrix cr A acts 
on the sublattice (A-B) space, while the matrix r 7 acts 
on the Fermi point "flavor" space, with A, 7 G {1,2,3}. 
Then 



(a x ,a y ) = (- 



and 



{P°,f3\(3 2 } = {r 



3 11 ^.1^21 



(6a) 



(6b) 



The disordered Dirac Hamiltonian Hjj given by equa- 
tion JSJ is the most general non-interacting 4-component 
Dirac Hamiltonian (up to irrelevant terms) consistent 
with the defining symmetries of the chiral sublattice 
model. The disordered lattice Hamiltonian defined in 
Eq. @ is invariant under time-reversal (T) and chiral (C) 
transformations. In the continuum theory, these trans- 
formations appear as 



T : * -> t 1 *, i 



(7) 
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FIG. 2: Hexagonal Brillouin zone appropriate to the triangu- 
lar sublattice. The independent zero energy Fermi points of 
the honeycomb model occur at ±fcjr; low-energy states have 
wavenumbers k such that |fc ^ &f < A <C &f ■ 

and 

C : * -> o- 3 (* f ) T , i -> -i (8) 

where T denotes matrix transpose. Note that both the 
time-reversal and chiral transformations are antiunitary 
in this language. We stress that in the 2 nd -quantized for- 
mulation discussed here, the "chiral symmetry" is rep- 
resented by the fact that the (second-quantized) Dirac 
Hamiltonian, given by Eq. JSJ), is invariant under the 
transformation defined by Eq. JSJ. This is easily seen to 
be equivalent to the first- quantized definition in Eq. (I}, 
above. — The Hamiltonian Ho in Eq. JSJ was originally 
derived in the context of the random bond TT-fhix model 
on the square lattice jiLiiLSi note that on the honeycomb 
lattice, we obtain the same effective field theory in zero 
magnetic flux. 

We take A^(x) and m l (x) to be Gaussian white-noise- 
distributed random variables, and we require that the 
ensemble-averaged system be invariant under honeycomb 
lattice rotations, translations, and reflections. These as- 
sumptions lead to the conditions 

A„{x) = m'(sB) = 0, (9) 



A^x)A v (x')=g A 2TrS^S(x-x'), (10) 

and 

m l (x) mi{x') = g2TrS ij S(x - x'). (11) 

In these equations the overbar denotes disorder averag- 
ing, and the factors of 2ir are conventional. The disorder- 
averaged theory is characterized by two parameters, the 
variances g and gA of the random mass and random vec- 
tor potential, respectively. 

An interpretation of the two types of disorder appear- 
ing in Eq. (0 in terms of the physics of the honeycomb 
lattice model is obtained by examining the structures of 
the matrices in Eq. ®, used to define them. The ran- 
dom masses couple the two Fermi nodes together in fla- 
vor space, while the random vector potential does not. 



Thus the random mass parameter g characterizes elastic 
backscattering between the two independent Fermi nodes 
(involving large crystal momentum transfers), while the 
random vector potential parameter gA characterizes in- 
tranode scattering involving small crystal momentum 
transfers. As discussed in Sec.H] g turns out to be related 
to the dimensionless DC resistance, which is independent 
of gA', on the other hand, gA dominates "one-particle" 
properties, including the density of states^ This separa- 
tion of one and two-particle properties is a characteristic 
feature of the chiral disorder classes. 



B. Adding interactions 

We complete the setup of our problem in this subsec- 
tion by adding to the disordered Dirac model [defined by 
Eq. generic short-ranged interactions that preserve 
the chiral (sublattice) symmetry. Consider the set of all 
four-fermion operators in the low-energy effective contin- 
uum theory that are purely local in space and time, and 
involve no derivatives. All such operators are irrelevant 
at the clean Dirac fixed point, but comprise the least ir- 
relevant set of local multiparticle interaction terms that 
we can write down. One then expects that these should 
be the most important local operators to study once we 
turn on the disorder. We can enumerate all independent 
four-fermion terms consistent with (i) the discrete T and 
C symmetries [Eqs. (JJJ and ©], and with (ii) the assump- 
tions of statistical invariance under honeycomb rotation, 
translation, and reflection transformations. It turns out 
that there are four such local operators; in this subsec- 
tion, however, we will only discuss two of them in detail. 
(The other two turn out to be even more irrelevant in the 
disordered phase than at the clean Dirac fixed point. In 
Sec. IIIII below, we will calculate the anomalous dimen- 
sions of all four operators.) 

Consider the local fermion density operators ua% — 
CAi°Ai and nBj — CgjCBj defined on the A and B sublat- 
tices of the honeycomb model. In the continuum Dirac 
theory, we define 

p(r) = -tf(r) (1 + a 3 ) ijj{r) ~ n A (r) (12a) 

and 

p(r) = ^(r) (1 - a 3 ) ^(r) ~ n B (r), (12b) 

where the identification with the sublattice operators 
holds up to terms that oscillate on the lattice scale. 

The two important four-fermion interactions may 
be written simply in terms of products of the long- 
wavelength sublattice density operators defined in 
Eq. (JT3J. They are pp(r) on one hand, and (p 2 + p 2 )(r) 
on the other. The operator pp, which involves a product 
of the long-wavelength densities from the two different 
sublattices, can be shown to result from the continuum 
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limit of, for example, nearest-neighbor density-density 
interactions 4 s . 

( n Ai ~ h)( n B] - |) ~ / d 2 rpp(r). (13) 

<«) 3 

The other important four-fermion local operator is 
p 2 + p 2 . Microscopically, the operators p 2 and p 2 are 
pure intrasublattice interactions, i.e. they do not cou- 
ple together particles on the two different sublattices. 
These are the only such four-fermion terms that we can 
write down. Note that invariance under y-axis reflections 
(which exchange the two sublattices — see Fig.0 requires 
p 2 and p 2 to appear in the combination p 2 + p 2 with 
identical coefficients. Thus p 2 + p 2 represents the most 
relevant piece of arbitrary short-ranged same-sublatticc 
(e.g. next-nearest-neighbor) interactions. 

In order to perform the average over disorder real- 
izations, we write a standard (2 + l)-dimensional (real- 



time-ordered (T-ordered) Grassmann path integral Zt 
to represent the disordered Dirac theory. We use the 
Schwinger-Keldysh method 37 to normalize this path in- 
tegral to unity in order to perform the disorder average; 
in this method, we exploit the fact that the normaliza- 
tion factor 1/Zt can be written simply as an anti-time- 
ordered (T-ordered) path integral Zf for the same theory. 
This construction works regardless of whether the model 
includes interparticle interactions or not. 

Specifically, we write the theory in terms of the 
Schwinger-Keldysh generating functional Z defined as 
(see e.g. Ref. Ha) 

Z = Z T Z T = J V^ a V^ a e l{So+Sl \ (14) 

where the non- interacting action Sd, containing the dis- 
order potentials, is given by 



Sp = ^ s a / dud 2 r ip a {u,r) [hw + is a f] sgn(w) + a p (id^ + A^f3°) - m l f3 l ] ip a (uj,r), (15) 

a=l,2 J 

while Si contains the four-fermion interactions 

Si = -J2 *a [ dtd 2 r[U(p 2 a + p 2 a )(t,r) + Wp aPa (t,r)+XOX(t,r) + YOUt,r)]. (16) 

a=l,2 ^ 

I 



The integration in Eq. (|14fl is over the Grassmann fields 
ip a and "tpa, where the ("Keldysh-") index a 6 {1,2} de- 
notes the T-ordered (a = 1) or T-ordered (a = 2) copy 
of the theory. The factor s a appearing in the action Sp 
(fT5)l is given by 

_ f 1, a = 1 (T-ordered), , _> 

Sa | -1, a = 2 (T-ordered). 

The term i s a r\ sgn(w) in Sd is the pole-prescription ap- 
propriate to T and T-ordered Green's functions, with 
X] — > + . In Eq. (JTSJ, we have set the Fermi velocity 
Vp = 1, and the various 4x4 matrices appearing in the 
theory were defined in Eq. © . Finally, the coupling con- 
stant h in Eq. (|15f) is a dynamic scale factor, and will be 
used to track the relationship between length and energy 
as we renormalize the theory. 

In Eq. I|16(l . we have introduced coupling constants U 
and W for the p 2 + p 2 and pp interaction operators, re- 
spectively. We denote by X and Y the coupling con- 
stants of a pair of additional interaction operators O x 
and O y that carry the same naive engineering dimension 
as p 2 + p 2 and pp. Although these additional operators 
are required to close the one-loop RG, we do not define 
them explicitly at this point (but see Sec. IIII F3 3|) . be- 



cause they will turn out to be strongly irrelevant in the 
disordered phase. To simplify notation, we define the 
vector of interaction couplings 

U ee {U, W, X, Y}. (18) 

The expressions for the Schwinger-Keldysh path inte- 
gral Z [Eq. and the action Sd + Si [Eqs. (JT5J and 
(|16fl ] constitute the completed formulation of the full, 
interacting and disordered fermion problem with chiral 
symmetry. Upon averaging over realizations of the disor- 
der potentials m l and A^, we find that the theory requires 
the specification of a total of seven coupling constants: 
the random mass and random vector potential variances 
g and gA, the four interaction strengths U, and finally 
the dynamic scale factor h. In the next section we per- 
form a perturbative one-loop RG analysis on this model. 
The principal new results are the scaling behaviors of the 
interaction couplings U as we look at lower and lower en- 
ergies, or, equivalently, larger and larger values of the 
flowing coupling gA- The reader interested only in our 
results may skip directly to Sec. IIVI where the scaling 
behaviors of the various couplings are simply stated and 
analyzed. 
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TABLE I: Continuum versions of honeycomb lattice model symmetry transformations in the L/R language. 



Symmetry 



Expression in terms of left(L)- and right(R)-moving components 



Chiral (C) 

Timc-rcvcrsal (T) 

x-Rcflcction" 

y-Rcflcction a 

Rotations' 1 

Translations" 

Electrical charge U(l) 



R(t) -> -^Li-t) 
R(t) -» -L(-t) 
R — > — L 



R 
R 
R 



"R 

">r 



L(t) -> -n x R{-t) 
L(t) -> -R(-t) 
L -> -fl 
L -» i/i 2 fl T 
L -» e _i(, L 
L -» e~ i() L 



fi(t) -> L(-t)^ 1 
R(t) -> £(-*) 

-> Re~ ia 
R -> He"* 9 * 3 



"Implicit transformation of the field argument r is implied. S 
So + Si denotes the action in Eg. 1141 . 



L(t) -> R(-t)^ 

£(t) -> fl(-i) 

£ -> fl 

£ -> -R T in 2 

I -» Le- i8 

£ -» £e i£) 

£ -> £e- ie " 3 



S -> S* 

s -> s* 



III. ONE-LOOP CALCULATION 

In this section we perform a one-loop renormaliza- 
tion group analysis on the interacting and disordered 
chiral fermion model defined in Sec. [H] We begin by 
re-writing the Grassmann Dirac field appearing in the 
Keldysh path integral [Eq. iflijl] in terms of "left(L)-" 
and "right (R)-moving" components, inspired by a sim- 
ilar decomposition made for the non-interacting system 
in Ref. 13. (This L/R decom pos ition provides a language 
similar to that used in Ref. [l/j to investigate the freez- 
ing transition in the dynamic critical exponent and the 
low-energy density of states.) We then perform a hard 
cutoff field theory renormalization group analysis to one 
loop. The L/R formulation allows for a particularly sim- 
ple and transparent treatment using diagrams. With a 
minimal effort we recover, as a check, the one-loop RG 
results for the disorder-only modelii*ii^4 [Eqs. 1391 
14 1 p ] . Finally, we calculate the anomalous dimensions of 
all allowed four-fermion interaction operators in the pres- 
ence of the disorder [Eqs. g3 HO ED EJ] . 



where the unitary matrix U is 

U = U p ^(l + ia 2 T*) ^(l-ia*). (20) 



Here U v 



(1/2) (1 



• ct • r ) is the permutation oper- 
ator, exchanging o a «-> r Q upon conjugation. This 
basis change induces a similarity transformation on 
the matrices defined in Eq. ijrjjl: the transformed ver- 
sions are a. = (a x ,a y ) — (o i ,o r ) and 
{/3°,/?\/? 2 }= {<rVV 3 r\r 2 }. 

Following Ref. ITU we define the L/R decomposition of 
the Dirac spinors: 



L2a 
R2a 

Rla 



(21) 



and 



A. Left (L) /Right (R) decomposition 



Ipa = S a [Rla Li a R 2a L 2a ] , 



(22) 



Consider the Keldysh action for the non-interacting 
disordered chiral system, given in Eq. (|15f) . We make the 
basis change [a = 1, 2 is the Keldysh-index; Eas. l|15ll6|) ] 



Ipa — » U 1p a , Ipa — > Ipa W , 



(19) 



where the first index i 6 {1,2} on the components R ia 
and Lia is a species index. Thus the fields R a ("right- 
movers" ) and L a ( "left-movers" ) are two-component ob- 
jects. In terms of the latter, the non-interacting action 
given by Eq. I|15|) may now be written very simply as 



S D = i Jdtd 2 r [h (Ra^dtLl] 



+ h {Ll^dtRa) + R a (d - A)R a + L a {5- A)L a 
+mR a L a + friL a Ra\ ■ 



(23) 



r 



where the Pauli matrix fi 1 is a matrix in the species space 
of the Fermion fields [which is the first index i of Eqs. H21|) 
and 1(22)) ]. 

In equation ll'.'il) we have suppressed the T and T- 
ordering pole prescriptions [Eq. (|15fl ] , and we have intro- 



duced the following complex notations: (i) d 



and d = d„ 



id 



(ii) A = Ay 

- im 1 and to 



iAr,. and A 



= 9 X 
A,. 



■Wy 

iA-t 

and (iii) m = to 2 

Turning to the interactions, we seek the set of all inde- 
pendent four-fermion operators local in space and time 
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with no derivatives; as discussed in subsection lll Bl these 
operators form a set of local perturbations least irrelevant 
at the clean Dirac fixed point. Any allowed interaction 
operator must be invariant under discrete chiral, time- 
reversal, and spatial reflection operations, and it must 
also transform as a singlet under continuous U(l) "elec- 
trical charge" , rotation, and honeycomb lattice transla- 
tion transformations 4^ We summarize the action of these 
symmetry transformations in the above L/R basis of low- 
energy fields in Table [I] Consistent with these require- 
ments, we find the four independent local four-fcrmion 
operators 



Ot EE 

and 
where 



- p 1 - p}^ - p 2 - p 2 

{Ra-^-Ra) (L a —L a ) + (R a —R a ) (L a —L a ), 



- p - p^ 

{Ra-^-Ra) {L a —L a ), 

(R a R a ) (L a L a ) 



— Ha ~ Fa, 



Pa 
Pa 



i{L a p 2 R T a ), 
-i{R T a p 2 L a ). 



(24a) 
(24b) 

(24c) 
(25) 



(26a) 
(26b) 



Equation l|26|l gives the L/R language expressions for the 
long-wavelength sublattice density operators defined by 
Eq. ■ The operator given by Eq. is one of the 
two four-fermion interaction terms that will turn out to 
be more relevant in the disordered critical phase than at 
the clean Dirac fixed point; the other, p a Pa, can be writ- 
ten in terms of the three operators defined in Eq. (|24H : 



Pa Pa 



-2 (<# + 0? + O c a 



(27) 



Note that the three operators in Eq. (|24|l . as well as the 
whole structure of the non-interacting theory [Eq. (|23|l ] 
except the "energy" terms (the terms involving time 
derivatives) are invariant under a U(l) transformation, 
which we refer to as a "sublattice" U(l), 



Ra > e * Rai L a 



L a , R a —> R a e % , L a — > L a e 



(28) 

(distinct from the "electrical charge" U(l) defined in 
Table I). This U(l) transformation is not a symmetry, 
because the energy terms in the non-interacting action 
carry non-zero (i.e. unit) charge. In fact, under this U (1), 
the two energy terms i(R a p 1 d t Lj) and i(L^p 1 dtR a ) 
carry equal and opposite unit charges, as do the lo- 
cal sublattice density operators p a and p a , while the 
four-fermion interaction operator [Eq. I|25|l ] is a sum 
of positive and negative doubly-charged same-sublattice 
components. 
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R 
L 



FIG. 3: Feynman rules for the diagrams in Fig. 2] 



In the non-interacting disordered fermion problem, 
moments of the local density of states (LDOS) operator 
are known to exhibit negative scaling dimensions in the 
disordered, critically delocalized phase^ii^ [The LDOS 
operator appears in the (2+l)D Keldysh formulation as 
the Dirac fermion bilinear ^> a {^>, T)i}) a ({jj, r), which is local 
in position space and frequency but not in time.] More 
precisely, it is the components of these LDOS moments 
with maximal charge under the sublattice U(l) trans- 
formation defined by Eq. (|28|l that carry the maximally 
negative scaling dimensions, i.e. these components are 
the most relevant perturbations to the zero-energy nearly 
conformal field theory. The functional RG treatment em- 
ployed in Rcf. 3 focuses just upon these most relevant 
components of the LDOS moments in order to extract 
the dynamic critical exponent in the regime above the 
freezing transition. In our interacting problem, we will 
find that the interaction operator = p\ + p\, which 
is maximally-charged under the "sublattice" U(l), is the 
most relevant of the investigated four-fcrmion interaction 
terms. Note that unlike the LDOS operator or the en- 
ergy terms in the non- interacting action (|23|l . the sum of 
sublattice density operators p a + p a cannot receive any 
corrections to scaling: this local spacetime operator is a 
component of the (2 + l)-dimensional conserved "elec- 
trical charge" U(l) current carried by the electrically- 
charged Dirac fermions [Table I]. 

Next we average the Keldysh functional Z [Eq. l|14Jl ] 
over realizations of the random vector potential A = A x + 
iA y and the complex random mass m, with A(x) A(x') = 
iirgA 5(x — x') and fh(x) m(x') — Aug 8{x — x'). The re- 
sulting non-interacting disorder-averaged action Sd may 
be written in its final form as a sum of two terms: 
Sd = Sr* + SSd, where 



S 



D 



and 



J dtd 2 r[R a L T a 



d hp}d t 




Ra 


hp}dt 8 




L a 



(29) 



SS D = J dtdtl d 2 r [g (R a L a )(t) {L b R b ){t') 

+g A (R a Ra)(t) (L b L b )(t')] (30) 

The average over disorder has coupled together fermion 
bilinears with different Keldysh indices at arbitrarily dis- 
tant times. 
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(a) 



(b) 





L Y L 




R R 



R 



R 



FIG. 4: Diagrams renormalizing the self-energy. 



B. Perturbative analysis 



(b) 

QM. 

O u : 



FIG. 5: Feynman rules for the diagrams in Figs. 15| The 
interaction vertex labeled O m denotes the sublattice U(l) 
charge- neutral operators defined by Eq. 1241 . while the ver- 
tex labeled O u depicts the p 2 component of the sublatticc- 
charged interaction operator p 2 + p 2 . 





We calculate one-loop corrections to the propagators in 
(|2*5)l and to the disorder vertices in l(3T)|> , and we compute 
the anomalous dimensions of the interaction operators 
given by Eqs. (P4"|) and (|25J). The dependence of these 
corrections on a hard momentum cutoff A will be used~ 
to derive the renormalization group equation summarized 
in Sect. IV. 



1. Self-energy 

In order to calculate the self-energy, it is useful to de- 
fine the fields \a and Xa via [Keldysh index a € {1,2}] 



and 



Xa 



Xa 



Ra 



[Ra L J a 



(31) 



(32) 



We can think of \a = {x-ya} y as a 2-component R/L- 
valued object; the index 7 6 {R,L} addresses the two 
components of the definition in Eq. I|31|) : each field R a 
or L a carries in addition a 2-component species index i, 
denned in Eq. (|2*T)l . (x-yia when all indices are displayed.) 

The propagator for the field Xa follows from Eq. (|29[1 . 
and is given by 



(Xa(uJ,k)Xa(u,k)} 



kk — h 2 (uj 2 + irjs a ) 



where k = k x - 



■ik y and k 



k.i 



k hojp: 1 
htu/j, 1 k 

(33) 

iky. In Eq. I)33|) we have 



restored the correct T and T-ordcring pole prescriptions 
in the two Keldysh species. We can also re- write the 
disorder-averaged action SSd given by Eq. (p?TTf) in terms 
of the x-ya fields ["T" denoting the transpose, and 7 e 
{L,R}]; 

M^=-Atti J dtdt'd 2 r [g {XRuXlaM {xlbXmW) 

-9 A (XRa XRa){t) (XLb XLb)(t')} ■ (34) 



We depict the disorder vertices g and gA with dashed 
lines, shown in Fig. U3 

Figure 01 shows the one-loop disorder contributions to 
the self-energy. In these diagrams the propagator (x 7 Xy) 
is represented as a directed solid line labeled with the 
endpoint indices (7, 7') € {R,L}. (Keldysh indices are 
suppressed.) We furnish the external x ne ld l e § s with 
the labels "R" and "L" to indicate which propagator is 
being corrected. 

The two diagrams appearing in Fig. ^a) corre- 
spond to the one-loop self-energy corrections to (xrXr)- 
Both of these diagrams vanish: the first diagram con- 
tributes nothing because it contains a free Keldysh index 
summation, 37 while the second diagram vanishes because 
the momentum integral of the purely left-moving propa- 
gator (xlXl) gives zero. Thus we find that there is no 
self-energy correction to the "R-R" sector of the x field 
propagator to one-loop order. 

By contrast, the diagrams contributing to the self- 
energy correction to (xrXl) shown in Fig.^Jb) lead to a 
non-zero renormalization of the dynamic scale factor h. 
The irreducible two-point vertex function in this channel 
to one loop is [cf. Eq. (|3^|) ] 



RL 



(oj,k)=ihojn 1 [l + 2(g + g A )lnA], (35) 



where A is an ultraviolet momentum cutoff and we did 
not write out ultraviolet-finite corrections. Differentiat- 
ing Eq. (|35|l with respect to In A and setting the result 
equal to zero leads to the expression^ 



d\nh 
dl 



2(5 + 9a), 



(36) 



where / = — In A is the logarithm of the RG length scale. 
The subscript "an" appearing on the left hand side of 
Eq. H36fl stands for "anomalous," as the above calculation 
has given us only the anomalous part of the beta function; 
in order to obtain the full one-loop flow equation for h, 
we must add the appropriate "engineering dimension" to 
Eq. (OU). This will be done in Sec.|IV| below. 

Although we have considered the propagator renormal- 
ization due only to the disorder, it can be easily seen 
that none of the local spacetime interaction operators 



FIG. 7: Diagrams renormalizing gA- 



FIG. 6: Diagrams renormalizing g. 

defined in Eqs. 1|24(1 and 125J1 can generate a non-zero 
one-loop contribution to the Dirac field self-energy; such 
a contribution would necessarily correspond to either the 
generation of a mass or a chemical potential shift, both 
forbidden by symmetry. 

2. Disorder self-renormalization 

Next we turn to the self-renormalization of the disor- 
der couplings g and gA- We will only need the diagonal 
components of the propagator given by Eq. I|33|) . so we 
revert back to the {R, L, R, L} component language [the 
species index i — 1,2, Eq. (|2*T1 will be suppressed]. 
The propagators we will need are 



(R a {w,k)R a (u;,k)) 



-ik 



kk — h 2 (uj 2 + irjSa) 



and 



(L a (ui,k)L a (oj,k)) 



—ik 



kk — h 2 (u> 2 + irjs a ) 



(37) 



(38) 



Note that both propagators are odd under the exchange 
k^-k. 

Now we will use directed solid and dotted lines to de- 
pict the (RR) and {LL) propagators, respectively. We 
will continue to use dashed lines to indicate the disorder 
vertices g and gA, as shown in Fig. [51 a). 

Figure shows the corrections to the g vertex in Eq. 
(|30|l . The In A contributions due to the first two diagrams 
exactly cancel those of the second two, leading to the 
anomalous part of the one-loop beta function 



0. 



(39) 



The cancelation is easily understood if one notes that 
all four diagrams in Fig. contain a closed momentum 
loop integration involving one of each diagonal propaga- 
tor given by Eqs. (|37|l and (|38|l . All four diagrams con- 
tribute the same absolute ultraviolet-divergent amplitude 
to the sum. Note that the loop momentum necessarily 
flows against one of the two directed chiral propagator ar- 
rows in the first pair of diagrams in Fig. El hi the second 
pair of diagrams, the loop momentum can be taken to 
flow in the same direction as both L- and R-propagator 
arrows. This simple structural difference produces a rel- 
ative minus sign, forbidding the coupling gA from feeding 



into g. The same mechanism will lead to cancelations in 
the autorenormalization of gA as well as in the computa- 
tion of the anomalous dimensions of some of the interac- 
tion operators. 

In Figure [7| we depict the one-loop disorder-disorder 
diagrams that renormalize the coupling gA- The first 
two diagrams, each proportional to g\, exactly cancel 
via the mechanism described above. The third diagram 
is non-zero, and leads to the irreducible four-point vertex 
function 



T gA = 5A + 2.g 2 lnA. 



(40) 



Differentiation with respect to / = — In A produces the 
anomalous part of the one-loop beta function 



d£A 

dl 



= 2g" 



(41) 



3. Interaction operator anomalous dimensions 

Finally, we turn to the one-loop calculation of the 
anomalous dimensions of the four interaction operators 
defined in Eqs. 1)24(1 and 125JI in the disordered (critical) 
phase characterized by finite g and gA- The results are 
listed in Eqs. 63 E3 113 EJ and Eqs. (GHH |6Qb| IfiUcl 
IflOdll . 

Consider first the three "sublattice-U(l)" charge- 
neutral interaction operators defined in Eq. I|24|l . We 
will use a wavy line to represent all three interaction 
operators simultaneously, as depicted by the vertex la- 
beled O m in Fig. Efb); these operators differ only in 
their species index structure. Figure [S] shows the one- 
loop renormalization of these interaction operators due 
to the disorder vertices g and gA- The labels i, j and fc, I 
appearing on the external legs denote the species compo- 
nents of the outgoing and incoming right and left-moving 
particles, respectively. The first four diagrams describe 
the (^-dressing of the interactions; these diagrams can- 
cel up to ultraviolet-finite contributions. The final two 
diagrams in Fig. [S] involve the g disorder vertex and add 
constructively. Let us denote the zeroth order irreducible 
four-point vertex function for the interaction operator 
O m via 



(42) 



In this equation, the vertex matrices M x depend upon 
the operator in question, e.g. the operator O a [Eq. (|24b 'l] 
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+ 




L-4 



+ 



2 
1 



2 
1 



+ 



FIG. 8: Diagrams renormalizing the interaction operators de- 
fined by Eq. fM . 



has two such vertex matrices, M 1 = p 1 /2 and M 2 = 
/i 2 /2. Incorporating the correction due to g, we find the 
one-loop irreducible four-point vertex function for the op- 
erator o M 



Tm irM = J2 M io M kl - 4 9 In A ]T M*M, 



,j [ ki- 



rn 



We have such a vertex function for each of the three 
operators defined in Eq. I|24[l . 

We rewrite Eq. 143|) using the Pauli matrix identities 



(44a) 



and 



(44b) 



(44c) 



where in this last identity a is summed over {1,2,3}. 
Using these identities with Eq. I|43|l shows us that the in- 
teraction operators mix under disorder-renormalization, 
leading to the matrix vertex equation 



(45) 



where the constant c = 2gln A. 

The eigenvalues of the matrix in Eq. I|45ll determine the 
anomalous dimensions of the RG eigenoperators, which 
are linear combinations of the interaction operators de- 
fined in Eq. i|24|l . These eigenvalues are 1 + Ag In A and 
1 — 4glnA; the latter is doubly degenerate. The eigen- 
operator corresponding to the eigenvalue 1 + 4g In A is 







1 2c 2c 




1 A 






c 1 — c c 




1 B 


r c 




c c 1 — c 




r° 



FIG. 9: Diagrams renormalizing the interaction operator p 2 



If we add this interaction operator to the action of the 
theory with a coupling constant W, as in Eq. I|16f) . then 
we find the anomalous part of the beta function for W 



dW 
~~dT 



= 4gW, 



(47) 



making the interaction coupling W more relevant in the 
(critical) disordered phase than at the clean Dirac fixed 
point. The other two eigenoperators share the RG eigen- 
value 1 — 4gln A; if we assign these operators the coupling 
constants X and Y, then we have 



and 



dX_ 



dY 



-igX, 



(48) 



(49) 



O a + O u + O c oc pp. 



(46) 



so that the couplings X and Y turn out to be more irrel- 
evant in the presence of disorder. 

Finally, we turn to the same-sublattice interaction op- 
erator O u given by Eq. I|25|l . We will use diagrams 
to compute the renormalization of the p 2 half of this 
operator; the other half p 2 must renormalize the same 
way. Using Eq. l(2rj|) . we write p 2 = — {L a p 2 R^) oc 
LiaLzaRiaRza- We will depict the interaction p 2 as a 
wavy line again joining right and left-mover lines, but 
now all such lines are directed outwards from the vertex — 
this half of the interaction operator is a source of "sub- 
lattice" U(l) current [see Fig. [Sfb)]. 

The diagrams describing the disorder-dressing of the 
p 2 operator are shown in Fig. [5] Numerical labels on the 
external legs again refer to the species index i, j = {1, 2}. 
The first four diagrams involve dressing by the i?^ vertex. 
The key point is that these four diagrams add construc- 
tively for this sublattice-charged interaction: gA feeds 
into the operator p 2 . The cancelation that occured for 
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the sublattice charge-neutral operators given by Eq. I|24|) 
does not occur here, because all four diagrams propor- 
tional to gA involve a loop momentum that necessarily 
flows against one of the two directed propagator arrows 
in that loop. The final two diagrams in Fig. [5] depict 
the g-renormalization of the interaction. Summing these 
contributions we obtain the irreducible four-point vertex 
function for the interaction operator O u = p 2 + p 2 to one 
loop 



l + 4(2.g A - 3 )lnA. 



(50) 



Adding O u to the action with coupling U as in Eq. (|16[) . 
we find the anomalous part of the beta function for the 
coupling U given by 



(51) 



IV. RESULTS AND DISCUSSION 

A. Scaling equations 

In this section we summarize the scaling behaviors 
of the coupling constants for the model defined by 
Eqs. O, {T3J), and (JTSJl in Sec. |TT| Seven coupling con- 
stants are needed to describe that model: these are the 
disorder variances g and gA, the interaction strengths 
U = {U, W, X, Y}, and the dynamical scale factor 
h. We may extract the full RG beta functions for these 
couplings by combining the results of Sec. lIIII with dimen- 
sional analysis, yielding the "engineering dimensions." 
Examine the action for the theory given by Eqs. (|15|) and 
(|16fl . We take the dependence of the energy "u>" upon 
the RG scale factor "Z" [the logarithm of the (spatial) 
length scale] to be 



dine 
~~dT 



z(l), 



(52) 



In Eq. l(52l . z is the (possibly scale-dependent) "dynamic 
critical exponent." We say that energy has the "engineer- 
ing dimension" 

[uj] = z (53) 

in inverse length units. We take the Fermi velocity vf = 
1 to be dimensionless; then, from Eq. l|15f) . the Dirac field 
has the engineering dimension 



W,r)] = 



1 + 2 



[i/j(u,r)] = 



(54) 



in two spatial dimensions. We then find that the disor- 
der strengths g and gA are dimensionless, the interaction 
couplings U share the engineering dimension 



[U] 



(55) 



and the dynamic scale factor h has the dimension 

[h] = 1 - z. (56) 

Combining the anomalous dimensions in Eqs. 1361 1 . 
Pty. gTJ, 1)47I49|) . and (EH from Sec.llTllwith the above 
engineering dimensional analysis, we find the one-loop 
flow equations 



e?ln h 
dl 

dg_ 

dl 



(l-z) + 2{g A +g) 
+ 0(g 2 ,gg A ,gU,gAU,U 2 ), (57) 
+ O(gU l9 3 ), (58) 



and 



dgA 
dl 



= 2g 2 + 0(g A U,g 3 ), (59) 



for the parameters that appear in the non-interacting 
sector of the theory, and 

^ = (-z + 8g A - 4g)U + O0 2 ,g 2 U, gg A U) (60a) 



^- = (- z + 4g)W + 0(U 2 ,g 2 W) 
^ r = (-z-4g)X + 0(U 2 ,g 2 X) 



as well as 



dY 
~~dl 



= (-z-4g)Y + 0(U 2 ,g 2 Y), 



(60b) 
(60c) 

(60d) 



for the interactions^ 

In Eqs. H57I59[> and (JHOJ), the RG function z(l) is so 
far undetermined. We may choose z any way we like; 
different choices correspond to different ways of rescaling 
energy in the (2+l)D field theory after performing the 
RG transformation. The natural choice 



dlnh 
dl 



0, 



(61) 



eliminates the parameter h from the theory. Enforcing 
Eq. (|61(l order by order in perturbation theory uniquely 
determines z; to one loop, we find that 



z(l) = l + 2(g A + g). 



(62) 



First, note (as a check) that we recover the known 
result o 11 ' 1 ?^ 4 (JEEl i 'EHl , and l(fi2")l for the non- interacting 
sector of the theory: the random mass variance g is 
purely marginal to one loop; by contrast, gA is driven to 
strong coupling by g ^ 0; the scale-dependent dynamic 
critical exponent z is fed by both g and gA- 

The scaling behaviors of the interaction couplings 
given in Eq. I|60|l are new and will be analyzed below. 
We stress that the terms maintained in the flow equations 
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(|58I60|) and l|62(l are only those required to determine the 
stability of the non-interacting phased 

Now we address an immediate concern: we are inter- 
ested in the disordered critical phase that occurs near 
zero energy as gA flows away from the clean Dirac fixed 
point (where g = gA = 0) to strong coupling; this flow 
always occurs for g ^ 0. The reader should therefore 
question the value of the purely perturbative one-loop 
results that we have presented above; naively, these low- 
est order results cannot be trusted to tell us anything 
about the strong coupling \aige-gA regime. 

Nevertheless, we can safely extend our results to large 
values of the random vector potential variance gA ^ 1; 
in order to do so, we rely upon several known exact re- 
sults for the disorder-only model 1315 . In the absence of 
interactions, the exact flow equation for g is strictly zero, 
while the exact flow equation for gA depends only upon 
g and vanishes for g = 0. The exact flow equation for 
the dynamical exponent z below the freezing transition 
is linear in gA, with non-linear ^-dependent coefficients. 
Equations l|59|) and l|62|) give the correct behaviors to low- 
est order in g. In this paper, we imagine tuning g <C 1; 
then our results in fact tell the whole story for the cou- 
plings g and gA, and for the dynamic critical exponent 
z. (Above the freezing transition, with gA > g\i we will 
need to be more careful — see below.) 



60b(l , we find using Eq. JSlJ 



U(l) 



"(0) 
"(0) 

Mi) 

and using Eq. (|4"7j) 



exp 



M>( dU 



exp[2(2 5 2 + (8.g^-4 5 )Z] (63) 



W(l) 



"(0) 
"(0) 



exp 



fdW 



. * v * 



exp [4g 1} 



(64) 



where we take g A -C 1 as the bare value of 17,4. In order 
to derive Eqs. 1|63[) and Q64II . we have integrated the flow 
equations for g and gA given by Eqs. and i|59|). 

To complete the story, we need the relationship be- 
tween the energy u>(l) and the RG parameter I] we ob- 
tain this by integrating Eq. I|62|l . We solve the resulting 
equation for / in terms of lnw, where u = w(0) is the 
bare energy scale. The result is 



(2 5 ) 5 



-l-2{g + g° A ) 



+ ^j[l + 2{g + g%)] 2 +8.g2 



In 



it 



(65) 



Indeed, 



the scaling dimensions of the important four- 

t2 1 „2 



fcrmion interaction operators p +p and pp, given to one 



loop in Eqs. H6Ua|) and (|60b|l . respectively (the terms lin- 
ear in U), must be exactly linear in gA, with higher-loop 
corrections only modifying their ^-dependences. These 
scaling dimensions combine the engineering dimension 
"— z" in these equations with the anomalous dimensions 
calculated in Sec. 1 1 i 1 1 The anomalous dimensions char- 
acterize the disorder-dressing of the interaction opera- 
tors at zero interaction coupling, and could alternatively 
be calculated in the (2 4- 0)D (non-interacting) version 
of the theory. Such a calculation in the (2 + 0)D the- 
ory can be performed non-perturbatively in gA^^; the 
anomalous dimension of a local operator in this (2 + 0)D 
theory is guaranteed to be linear in gA with a coefficient 
determined by the "sublattice" [Eq. JUJl] and "transla- 
tion" [see TableU] U(l) charges carried by that operator. 
Therefore our results in Eqs. (|60a|l and l)60b|l are in fact 
valid to all orders in gA- Eqs. (|6Ua[) and (|60b(l incorpo- 
rate, in addition, the lowest order corrections in g. 

In the following, we will analyze our results. We note 
from Eq. (|SU|l that the interaction couplings U and W are 
enhanced by the disorder, while X and Y are (weakly) 
suppressed relative to the clean Dirac fixed point. We 
therefore focus our attention upon the former pair of 
interactions. Integrating the flow equations (|60a|l and 



where Q = u>(l) > u) is an arbitrary finite reference 
energy. We will state all of our results in terms of the bare 
energy u> — > 0, so it is useful to note that the coupling gA 
evolves according to 



g A (l) ~ 2g 2 l + gZ, 



(66) 



and we use Eq. (|65|l to relate I to u>. 

Now we examine the scaling behaviors of the interac- 
tion strengths given by Eqs. I|63|) and (|64|) in the follow- 
ing three dynamical regimes, which we distinguish by the 
strength of the vector potential disorder 17^4: (1) near the 
clean Dirac fixed point (gA ~ g% <C 1), (2) at interme- 
diate disorder below the freezing transition (gA — > gX), 
and finally (3) above the transition at asymptotically 
strong disorder (gA 3> gX)- The freezing transition 
occursi 5 - at gA = 1 + 0(g) for small g. The idea is to 
tune both g and gA small at some arbitrary finite energy 
scale iv = fi; this places us close to the clean Dirac fixed 
point, where we know the interactions characterized by 
U and W are both strongly irrelevant. We then examine 
successively lower and lower energies lu — » 0; as we do 
so, the RG sends gA to larger and larger values, and we 
expect scaling behavior of the interaction operators to 
change as we approach the delocalized, critical phase at 
zero energy. 

We now discuss in turn the regimes (1), (2), and (3) 
defined above: 
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1. Near the clean Dirac fixed point: < 1 



We begin near the clean Dirac fixed point with 
g 2 |ln(u;/Q)| <C 1. We can then approximate the dynamic 
scaling relationship l|65[) as 



LO 



(1 - 2.g - 2g° A ) , 



(67) 



which corresponds to a dynamic critical exponent z ~ 1 + 
2g+2g A ] , implying a vanishing density of states. Inserting 
Eq. 1|67|) into Eq. I|66[) gives the consistency check gA <C 1. 

Using Eq. (|cJ7|) , the scaling equations l|rj3*|) and (|rJl|) for 
the couplings U and W become 



(68) 



U/U° 



and 




W 



(- 



l-4g 



(69) 



As expected, near the clean Dirac fixed point both inter- 
actions are strongly irrelevant as we look at decreasing 
energy scales. 



2. "Intermediate" disorder: gA — > 9 a 

For larger values of the disorder gA, or equivalently, 
smaller energies with g 2 |ln(u;/Q)| ~ 1, Eq. (|65|l predicts 
a crossover in the dynamic scaling behavior from Eq. i|67|) 
to 



I 



1 

2fJ 



In- 



n 



const., 



leading to the scaling behaviors 

3 

exp 



LO 



9 



lag 



corresponding to the local interaction p 2 



and 



LO 











bag 




2 V 2 


f2 





(70) 



(71) 



(72) 



corresponding to the local interaction pp. Equation l|71|) 
is one of the primary results of our paper. 

Equation (|71|l tells us that the same-sublattice four- 
fermion interaction p 2 + p 2 has now become extremely rel- 
evant towards the low-energy limit lo — > 0, while the local 
interaction pp that couples different sublattices together 
remains irrelevant, albeit enhanced when compared to 
the clean Dirac point. Our result means that the low- 
energy, delocalized phase of 2D non-interacting fcrmions 
made possible by the very special chiral sublattice sym- 
metry is unstable to same-sublattice (e.g. next-nearest- 
neighbor) interaction effects. On the honeycomb lattice, 



FIG. 10: Scaling evolution of U(gA) for several fixed values of 
g. This figure tracks the irrelevance/relevance of the coupling 
U associated with the same-sublattice interaction p 2 +p 2 as we 
renormalize to stronger vector potential disorder gA — * oo, 
or equivalently, smaller energies u) — > 0. U is normalized 
relative to its the bare value U° at reference energy scale S7, 
and we have tuned g\ = 0. The non-interacting, critically 
delocalized phase resides at u> = with gA = oo. Initially 
irrelevant near the clean Dirac fixed point (gA 1), U 

ultimately flows off to strong coupling at larger gA, taking the 
model away from the non-interacting, zero energy delocalized 
phase. 



next-nearest-neighbor interactions tend to favor charge- 
ordering on one triangular sublattice for U < 0, while 
for U > these interactions are frustrated by the non- 
bipartite nature of the triangular sublattice. 

Figure ^| shows the RG evolution of the same- 
sublattice interaction strength U with increasing length 
scale as a function of the running random vector poten- 
tial variance gA for several different values of the fixed 
random mass variance g [see Eq. I|73|l below]. This plot 
describes the following situation: we turn on an infinitesi- 
mally small interaction strength U and set the bare vec- 
tor potential disorder g A — at energy u> = Cl. For small 
g this places us near the clean Dirac fixed point. Then 
we run the RG, flowing to smaller energies lo — > and 
stronger disorder gA — * oo, and we observe the evolution 
of the coupling U. When the dynamic scaling behavior 
crosses over from Eq. I|67f) to Eq. 1)71) [I. this interaction 
strength begins to diverge strongly. The analytical ex- 
pression plotted in Fig. 1101 is given by 



U(g A ) ~ exp 



3 

2? 



1 



9A 



(73) 



Equation l|73() may be obtained by combining Eqs. I|63|l 
iJHSJl, and (|ofj|l . setting g A = 0. The value of gA 
where U becomes relevant varies with g: for smaller val- 
ues of g, the crossover occurs at smaller values of g A 



9% 
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From Eq. J73J), we have 



9l 



6 



a- 



(74) 



Note, however, that the transition at small g occurs at 
much lower energies [since we need g 2 |ln(a;/f2)| to grow 
from zero {uj — fi) to order one {uj ~ Slexp(— 1/g 2 ) <C 
0)]. The freezing transition in the dynamic scaling be- 
havior is supposed to occur— at g% ~ 1 + 0{g) for small 
g; Eq. (|74|l suggests that for g <C 1, the interaction {J 
becomes relevant for we/Z 6eZoio g^, potentially pre- 
venting the RG from ever reaching the non-interacting 
freezing transition. 

The instability to interaction effects that we observe 
here is perhaps not so surprising when one considers 
the strongly-divergent form of the low-energy density of 
states in this "intermediate" disorder regime given by 
Eq. (J2J. The exponent is a = 1/2 in this equation. Wc 
have a picture of a large number of exotic, multifrac- 
tal states crammed into a very narrow energy window; 
we might therefore expect that turning on interparticle 
interactions should produce a very strong effect. The 
obvious question, then, is why does the same instability 
not occur for weak interactions in the corresponding ID 
problem? 

The non-interacting disordered fermion problem with 
chiral symmetry [Eq.Q] in one dimension may be for- 
mulated as random spin 1 /2 quantum XY chain— This 
random XY chain is asymptotically solvable by a strong 
randomness RG procedure^ 6 ^ wherein the spins con- 
nected by the strongest random bonds are systemati- 
cally decimated out of the chain, replaced by ever weaker 
bonds connecting the surviving spins. The infrared limit 
of this RG is the so-called random singlet phase 16 , which 
is similar—, in some respects, to the non-interacting ver- 
sion of the 2D model considered in this paper. The ad- 
dition of weak, homogeneous, short-ranged Ising SfS? 
couplings to the random XY spin chain, corresponding 
to density-density interactions in the fermion model, does 
not destabilize this disordered phase; in the strong ran- 
domness RG, these additional Ising bonds become ever 
weaker as spins coupled with strong XY bonds are re- 
moved from the chain. Although some progress^ has 
been made in applying these ideas to the non-interacting 
2D chiral Dirac theory, we do not know whether this 
methodology may shed further light on the instability to 
same-sublattice interaction effects that we have found in 
this paper. Presumably, the higher degree of connected- 
ness in the 2D model changes the way that the interac- 
tions are renormalized relative to the ID case. 

Recall that Eq. i|71|) only gives the scaling behavior of 
the interaction p 2 + p 2 (in the intermediately-disordered 
regime); in Eqs. (|58I6UII . we have not written explicitly 
the one-loop renormalization effects of U upon the dis- 
orders strengths g and gA or the autorenormalization ef- 
fects of U upon itself. Although we will omit the details 
here, we have performed the complete one-loop RG cal- 
culation and found that all of these corrections vanish 



to one loop. (This eliminates the possibility of finding a 
new, interacting critical fixed point to lowest order in U.) 
We also know that U cannot renormalize h at one loop, 
as discussed in the paragraph below Eq. H36|) in Sec. II I II 
and thus cannot influence the dynamic scaling behavior 
given by Eq. l|65l) . 

The scaling equations (|71|l and l|72|l derived above ap- 
ply equally to the sublattice hopping model on the square 
lattice with 7r-flux model per plaquette and real random 
bond disorder— (The random Hamiltonian is again in- 
variant under both chiral and time-reversal symmetries.) 
The operators p and p have the same interpretations as 
local sublattice density operators because their defini- 
tion [Eq. Q12p] relies only upon the definition of the chi- 
ral symmetry in the continuum Dirac theory, given by 
Eq. JSJ). The sublattices of the bipartite square lattice are 
themselves bipartite square lattices, and next-nearest- 
neighbor interactions tend to produce charge-ordering on 
one sublattice for U < 0, and charge stripes for U > 0. 



3. Strong disorder: gA ^> g% 

We can imagine that the RG flow reaches successfully 
the regime above the freezing transition gA > g% before 
U grows large enough to invalidate our scaling analysis, 
say by tuning the initial interaction strength very small. 
The ultimate dynamic scaling behavior at asymptotically 
low energies is different from that of Eq. I|65l) : instead, for 
9 A > 9ai the relationship between the RG scale I — > oo 
and the bare energy scale ui — ► is given b y 14 ' 15 



1 

<9) 



In- 



(75) 



with c(g) oc g + 0{g 2 ). Using the same logic as above, 
this leads to the scaling behaviors 



and 



C/^-exp 



1 

<9) 



1 

<9) 



n 



(76) 



(77) 



Equation (|76() is another primary result of our paper. 

The main point that we want to stress is that Eq. I|76|) 
shows that U is ultimately even more relevant beyond 
the freezing transition, so that regardless of whether RG 
takes us into the regime above the freezing transition be- 
fore U diverges, we always find that the non-interacting 
description is unstable to the same-sublattice interaction 
effects characterized by the operator p 2 + p 2 . We are able 
to extend our analysis to this large gA regime because the 
key results leading to Eqs. (|76|l and (|77|l can be obtained 
from the (2 + 0)D theory non-perturbatively in 17^4: as 
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discussed above Eqs. (|63|l and Q64J1. the 1-loop flow equa- 
tions for the disorder couplings g and [Eqs. (|5*S|l and 
(jSnil] j and for the interaction scaling dimensions [Eq. ffiH]!] 
contain the exact (^-dependence (to all orders in gA), 
and we combine these results with the non-perturbative 
dynamic scaling relation beyond the freezing transition 
[Eq. (O, taken from Ref. [jjj. This leads to Eqs. fflft 
and J77J- 

B. Summary and outlook 

Within the context of our perturbative one-loop renor- 
malization group calculation, we have found that same- 
sublattice interactions become strongly enhanced as one 
examines the TRI random chiral symmetric Dirac model 
at ever lower energy scales, or equivalently, larger values 
of the running disorder coupling gA . We emphasized that 
while the calculation presented in this paper appears per- 
turbative in both of the disorder strengths g and gA, we 
can extend our results to arbitrarily large values of gA 
because of the knowniii^ non-perturbative structure of 
the non-interacting (nearly conformal) disordered Dirac 
theory. 

The validity of our stability analysis is restricted to 
g -C 1; because g does not evolve under the RG in the 
non-interacting phase, we can tune g as small as we wish. 

The same mechanism responsible for the instability to 
interaction effects drives the divergence of the density of 
states, and ultimately the freezing transition, in the non- 
interacting model. We have shown that this instability 
occurs regardless of whether we reach the regime beyond 
the freezing transition or not. 

To lowest order in the interaction strength, we see no 
signs of the possibility of a new, interacting critical fixed 
point. We therefore anticipate that the strongly rel- 
evant same-sublattice interactions drive the disordered 
honeycomb model toward some kind of Mott insulating 



phase. We pointed out that the same picture obtains 
for the square lattice 7r-fiux modelii. In the introduc- 
tion we dubbed the observed evolution of the clean but 
interacting Dirac fixed point, weakly perturbed by chiral- 
symmetric disorder g at some finite energy scale to = Q, 
toward a Mott insulating ground state in the limit of zero 
energy u> — > 0, a "disorder-driven Mott transition." In 
the limit of small g, this instability occurs only at very 
small energies; the scale at which the same-sublattice 
interactions become relevant occurs below the freezing 
transition for small g, and is roughly the same scale at 
which the density of states begins to diverge, as in Eq. @ 
with a = 1/2. 

An obvious extension of the present work would be to 
try to approach the interaction effects beyond the regime 
linear in the strengths of the interactions; indeed, the 
calculation presented in this paper treats the interac- 
tions only to lowest order, addressing only the question 
of the stability of the non-interacting description to the 
inclusion of interactions. The natural framework for such 
an extended investigation is the generalized Finkel'stein 
NLaM 17 , which grafts many-particle interactions from 
the disordered fermion theory directly into Wegner's low 
energy effective theory of interacting diffusion modes. An 
advantage of the Finkel'stein NLctM is that the pertur- 
bative RG is an expansion only in the dimensionless DC 
resistance; therefore, at least formally, even a one-loop 
calculation treats the interaction couplings to all orders. 
A calculation along these lines is currently in progress^. 
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